#include <opencv2\core.hpp>
#include <opencv2\ml.hpp> 
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <sstream>
#include <iterator>  
#include<cstdlib>
using namespace std;
using namespace cv;
using namespace cv::ml;

void creatMat(Mat &x, Mat &y, String fileName) {
	int line_count = 0;
	char buffer[256];
	ifstream in(fileName);
	if (!in.is_open()) {
		cout << "Error opening file"; exit(1);
	}
	while (!in.eof())
	{
		in.getline(buffer, 100);
		stringstream stream;
		stream << buffer;
		string temp_s;
		stream >> temp_s;
		double num1 = atof(temp_s.c_str());
		if (num1 == 1.0) {
			y.at<double>(0, line_count) = num1;
		}
		while (stream >> temp_s) {
			int index = temp_s.find(':');
			string temp1_s = temp_s.substr(0, index);
			double t1 = atof(temp1_s.c_str());
			string temp2_s = temp_s.substr(index + 1, temp_s.length());
			double t2 = atof(temp2_s.c_str());
			//cout << t1 << "::"<<t2<<" ";

			x.at<double>(t1 - 1, line_count) = t2;
		}
		line_count++;
	}
}


void initial_parermaters(Mat &w, double &b, int n1, int n0) {
	w = Mat::zeros(n1, n0, CV_64FC1);
	b = 0.0;

	//double temp = 2 / (sqrt(n1));
	double temp = sqrt(6.0 / (double)(n1 + n0));
	RNG rng;
	for (int i = 0; i < w.rows; i++) {
		for (int j = 0; j < w.cols; j++) {
			w.at<double>(i, j) = rng.uniform(-temp, temp);//xavier初始化
														  //w.at<double>(i, j) = 0;
		}
	}
}

void sigmoid(const Mat &original, Mat &response) {
	response = original.clone();//防止维度不同
	double temp;
	for (int i = 0; i < original.rows; i++) {
		for (int j = 0; j < original.cols; j++) {
			temp = original.at<double>(i, j);
			response.at<double>(i, j) = 1.0 / (1.0 + exp(-temp));
		}
	}
}

void relu(const Mat &original, Mat &response) {
	response = original.clone();//防止维度不同
	for (int i = 0; i < original.rows; i++) {
		for (int j = 0; j < original.cols; j++) {
			if (original.at<double>(i, j) < 0) {
				response.at<double>(i, j) = 0.0;
			}
		}
	}
}

void change_log(const Mat &original, Mat &response) {
	response = original.clone();//防止维度不同
	double temp;
	for (int i = 0; i < original.rows; i++) {
		for (int j = 0; j < original.cols; j++) {
			temp = original.at<double>(i, j);
			response.at<double>(i, j) = log(temp);
		}
	}
}

void linear_activation_forward(Mat &a_prev, Mat &a, Mat &w, double &b, string activation) {
	cv::Mat z;
	if (activation == "sigmoid") {
		z = (w*a_prev) + b;
		//cout << w.rows<<","<<w.cols<<" " << a_prev.rows<<","<<a_prev.cols<<endl;

		sigmoid(z, a);
	}
	else if (activation == "relu") {
		z = (w*a_prev) + b;
		//cout << w.rows << "," << w.cols << " " << a_prev.rows << "," << a_prev.cols << endl;

		relu(z, a);
	}
}


float compute_cost(const Mat &y, const Mat &a) {
	double cost = 0.0;
	cv::Mat temp1 = cv::Mat::zeros(a.rows, a.cols, CV_64FC1);
	cv::Mat temp2 = cv::Mat::zeros(a.rows, a.cols, CV_64FC1);
	change_log(a, temp1);
	change_log(1 - a, temp2);
	cv::Mat loss;
	loss = y.mul(temp1) + (1 - y).mul(temp2);
	cost = (-1.0 / y.cols)*sum(loss)[0];
	return cost;
}


void activation_backward(const Mat &a, const Mat &da, Mat &dz, string activation) {
	if (activation == "sigmoid") {

		dz = da.mul(a.mul(1 - a));
	}
	else if (activation == "relu") {
		dz = da.clone();//保证维度相同
		for (int i = 0; i < a.rows; i++) {
			for (int j = 0; j < a.cols; j++) {
				if (a.at<double>(i, j) <= 0) {
					dz.at<double>(i, j) = 0.0;
				}
			}
		}
	}

}

void linear_backward(const Mat &da, const Mat &a, const Mat &a_prev, Mat &w, double &b, Mat &dw, double &db, Mat &da_prev, const int m, const double learning_rate, string activation) {
	cv::Mat dz;
	activation_backward(a, da, dz, activation);//激活函数的反向传播

	dw = (1.0 / m)*dz*a_prev.t();
	db = (1.0 / m)*sum(dz)[0];
	da_prev = w.t()*dz;
	w = w - (learning_rate * dw);
	b = b - (learning_rate * db);
}

//计算分类精度
double calculateAccuracyPercent(const Mat &original, const Mat &predicted)
{
	return 100 * (float)countNonZero(original == predicted) / predicted.cols;
}
int main()
{
	//读取数据，并保存到矩阵中
	cv::Mat train_x = cv::Mat::zeros(123, 1605, CV_64FC1); // 全零矩阵
	cv::Mat train_y = cv::Mat::zeros(1, 1605, CV_64FC1); // 全零矩阵
	string filename1 = "C:\\Users\\thinkpad\\Desktop\\data\\a.txt";
	creatMat(train_x, train_y, filename1);
	cv::Mat test_x = cv::Mat::zeros(123, 30956, CV_64FC1); // 全零矩阵
	cv::Mat test_y = cv::Mat::zeros(1, 30956, CV_64FC1); // 全零矩阵
	string filename2 = "C:\\Users\\thinkpad\\Desktop\\data\\a-test.txt";
	creatMat(test_x, test_y, filename2);

	int layer_dims[] = { 123,32,8,1 };//去掉首，3层神经网络
	cv::Mat w1, w2, w3 ;
	double b1, b2, b3 ;
	initial_parermaters(w1, b1, layer_dims[1], layer_dims[0]);
	initial_parermaters(w2, b2, layer_dims[2], layer_dims[1]);
	initial_parermaters(w3, b3, layer_dims[3], layer_dims[2]);

	cv::Mat a1, a2, a3;
	cv::Mat da0, da1, da2, da3;


	cv::Mat dw1, dw2, dw3;
	double db1, db2, db3;

	int num_iterations = 8000;//迭代次数
	double learning_rate = 0.1;//学习速率
	int m = train_y.cols;

	//线性传播
	ofstream out("C:\\Users\\thinkpad\\Desktop\\out.txt");
	
	for (int i = 1; i <= num_iterations; i++) {
		linear_activation_forward(train_x, a1, w1, b1, "relu");
		linear_activation_forward(a1, a2, w2, b2, "relu");
		linear_activation_forward(a2, a3, w3, b3, "sigmoid");	
		
		da3 = (1-train_y)/(1-a3)-(train_y/a3);//loss对a的导数														  //cout << da5;
	
		linear_backward(da3, a3, a2, w3, b3, dw3, db3, da2, m, learning_rate, "sigmoid");
		linear_backward(da2, a2, a1, w2, b2, dw2, db2, da1, m, learning_rate, "relu");
		linear_backward(da1, a1, train_x, w1, b1, dw1, db1, da0, m, learning_rate, "relu");
		
		double cost = compute_cost(train_y, a3);
		if (out.is_open())
		{
			out << i << ":" << cost<<endl;
			
		}
		if (i % 100 == 0) {
			cout << "After " << i << "*100 iterations cost" " : " << cost << endl;

		}
	}
	train_y.convertTo(train_y, CV_32S);  //转换为整型
	a3.convertTo(a3, CV_32S);  //转换为整型
	cout << "train accuracy: " << calculateAccuracyPercent(train_y, a3) << "%" << endl;

	linear_activation_forward(test_x, a1, w1, b1, "relu");
	linear_activation_forward(a1, a2, w2, b2, "relu");
	linear_activation_forward(a2, a3, w3, b3, "sigmoid");

	test_y.convertTo(test_y, CV_32S);  //转换为整型
	a3.convertTo(a3, CV_32S);  //转换为整型
	cout << "test accuracy: " << calculateAccuracyPercent(test_y, a3) << "%" << endl;
	system("pause");
}